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SUMMARY 

We derive electromagnetic finite elements based on a variational principle that uses the electro- 
magnetic four-potential as primary variable. This choice is used to construct elements suitable for 
downstream coupling with mechanical and thermal finite elements for the analysis of electromag- 
netic/mechanical systems that involve superconductors. The key advantages of the four-potential 
are: the number of degrees of freedom per node remain modest as the problem dimensional- 
ity increases, jump discontinuities on interfaces are naturally accomodated, and static as well 
as dynamics are included without any a priori approximations. The new elements are tested 
on an axisymmetric problem under steady-state forcing conditions. The results are in excellent 
agreement with analytical solutions. 
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1. INTRODUCTION 


The present work is part of a research program for the numerical simulation of electro- 
magnetic/mechanical systems that involve superconductors. The simulation involves the 
interaction of the following four components: 

(1) Mechanical Fields : displacements, stresses, strains and mechanical forces. 

(2) Thermal Fields : temperature and heat fluxes. 

(3) Electromagnetic (EM) Fields : electric and magnetic field strengths and fluxes, cur- 
rents and charges. 

(4) Coupling Fields : the fundamental coupling effect is the constitutive behavior of 
the materials involved. Particularly important are the metallurgical phase change 
phenomena triggered by thermal, mechanical and EM fields. 

1.1 Finite Element Treatment 

The first three fields (mechanical, thermal and electromagnetic) are treated by the finite 
element method. This treatment produces the spatial discretization of the continuum into 
mechanical, thermal and electromagnetic meshes of finite number of degrees of freedom. 
The finite element discretization may be developed in two ways: 

(1) Simultaneous Treatment. The whole problem is treated as an indivisible whole. The 
three meshes noted above become tightly coupled, with common nodes and elements. 

(2) Staged Treatment. The mechanical, thermal and electromagnetic components of the 
problem are treated separately. Finite element meshes for these components may 
be developed separately. Coupling effects are viewed as information that has to be 
transferred between these three meshes. 

The present research follows the staged treatment. More specifically, we develop finite 
element models for the fields in isolation, and then treat coupling effects as interaction 
forces between these models. This “divide and conquer” strategy is ingrained in the parti- 
tioned treatment of coupled problems [4,16], which offers significant advantages in terms of 
computational efficiency and software modularity. Another advantage relates to the way 
research into complex problems can be made more productive. It centers on the obser- 
vation that some aspects of the problem are either better understood or less physically 
relevant than others. These aspects may be then temporarily left alone while efforts are 
concentrated on the less developed and/or more physically important aspects. The staged 
treatment is better suited to this approach. 

1.2 Mechanical Elements 

Mechanical elements for this research have been derived using general variational principles 
that decouple the element boundary from the interior thus providing efficient ways to work 
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out coupling with non-mechanical fields. The point of departure was previous research into 
the free-formulation variational principles reported in Ref. [5]. A more general formulation 
' for the mechanical elements, which includes the assumed natural strain formulation, was 
established and reported in Refs. [5,6,14,15]. New representations of thermal fields have 
not been addressed as standard formulations are considered adequate for the coupled-field 
phases of this research. 

2. ELECTROMAGNETIC ELEMENTS 

The development of electromagnetic (EM) finite elements has not received to date the 
same degree of attention given to mechanical and thermal elements. Part of the reason 
is the widespread use of analytical and semianalytical methods in electrical engineering. 
These methods have been highly refined for specialized but important problems such as 
circuits and waveguides. Thus the advantages of finite elements in terms of generality have 
not been enough to counterweight established techniques. Much of the EM finite element 
work to date has been done in England and is well described in the surveys by Davies [l] 
and Trowbridge [21]. The general impression conveyed by these surveys is one of an un- 
settled subject, reminiscent of the early period (1960-1970) of finite elements in structural 
mechanics. A great number of formulations that combine flux, intensity, and scalar po- 
tentials are described with the recommended choice varying according to the application, 
medium involved (polarizable, dielectric, semiconductors, etc.) number of space dimen- 
sions, time-dependent characteristics (static, quasi-static, harmonic or transient) as well 
as other factors of lesser importance. The possibility of a general variational formulation 
has not apparently been recognized. 

In the present work, the derivation of electromagnetic (EM) elements is based on a vari- 
ational formulation that uses the four-potential as primary variable. The electric field is 
represented by a scalar potential and the magnetic field by a vector potential. The for- 
mulation of these variational principle proceeds along lines previously developed for the 
acoustic fluid problem [7,8]. 

The main advantages of using potentials as primary variables as opposed to the more 
conventional EM finite elements based on intensity and/or flux fields are, in order or 
importance: 

1. Interface discontinuities are automatically taken care of without any special interven- 
tion. 

2. No approximations are invoked a priori since the general Maxwell equations are used. 

3. The number of degrees of freedom per finite element node is kept modest as the 
problem dimensionality increases. 

4. Coupling with the mechanical and thermal fields, which involves derived fields, can 
be naturally evaluated at the Gauss points at which derivatives of the potentials are 
evaluated. 

Following a recapitulation of the basic field equations, the variational principle is stated. 
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The discretization of these principle into finite element equations produces semidiscrete 
dynamical equations, which are specialized to the axisymmetric case. These equations are 
validated in a simulation of a cylindrical conductor wire. 


3. ELECTROMAGNETIC FIELD EQUATIONS 


3.1 The Maxwell Equations 

The original Maxwell equations (1873) involve four spatial fields: B, D, E and H. Vec- 
tors E and H represents the electric and magnetic field strengths (also called intensities), 
respectively, whereas D and B represent the electric and magnetic flux densities, respec- 
tively. All of these are three-vector quantities, that is, vector fields in three-dimensional 
space (xx = x, X 2 = y, X3 = z): 


E = l 


E 2 
[e 3 


D = 




( 1 ) 


Other quantities are the electric current 3-vector j and the electric charge density p (a 
scalar). Units for these and other quantities of interest in this work are summarized in 
Tables 1-2. 

With this notation, and using superposed dots to denote differentiation with respect to 
time t, we can state Maxwell equations as* 


(2) 


The first and second equation are also known as Faraday’s and Ampere-Maxwell laws, 
respectively. 

The system (2) supplies a total of eight partial differential equations, which as stated are 
independent of the properties of the underlying medium. 

S.2 Constitutive Equations 

The field intensities E and H and the corresponding flux densities D and B are not 
independent but are connected by the electromagnetic constitutive equations. For an 
electromagnetically isotropic, non-polarized material the equations are 


B + V x E = 0, 

V xH-D = j, 

II 

P 

> 

V • B = 0 . 


(3) 


* Some authors, for example Eyges [2], include 4 n factors and the speed of light c in the Maxwell 
equations. Other textbooks, e.g. [19,20], follow Heaviside’s advice in using technical units 
that eliminate such confusing factors. 


B = /xH, D = cE 
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Tkble 1 Electric and Magnetic Quantities 


Quantifies 

Symbol MKS- Weber Units 

Electric chavge density 

9 

coulomb/m 3 

Electric field intensity 

E 

newton/ coulomb 

Electric flux density 

D 

coulomb/m 3 

Electric resistance 

R 

ohm 

Electric conductivity 

? 

mho 

Displacement current density 

D 

coulomb/ (sec.m 3 ) 

Susceptibility" 

€ 

coulomb/ (joule.m) 

Current 

j 

coulomb/sec 

Magnetic field intensity 

H 

newton/weber or amperes/m 

Magnetic fiuoc density 

B 

weber /m 3 

Magnetic permeability! 

M 

weber/ (joule un) or henry /m 

* Also called capacitivity and permittivity 
f Also called inductivity 


Table 2 Equivalences Between Various MKS-Glorgi Units 


1 newton = 1 kg.m/sec 3 
1 joule = 1 newt on. m 
1 watt = 1 joule/sec 
1 coulomb = 1 ampere.sec 
1 ohm = 1 volt/ampere 
1 farad = 1 coulomb/ volt 
1 henry = 1 (voltjec) /ampere 
1 weber = 1 volt.sec 
1 mho = 1 ohm -1 


where n and e are the permeability and susceptibility, respectively, of the material!. These 
coefficients are functions of position but (for static or harmonic fields) do not depend on 
time. In the general case of a non-isotropic material both n and e become tensors. Even 
in isotropic media \i in general is a complicated function of H; in ferromagnetic materials 
it depends on the previous history (hysteresis effect). 

In free space n = no and c = eo, which are connected by 


cl = 


MO^O 


( 4 ) 


where Co is the speed of fight in a free vacuum. In MKS-A units, c 0 = 3.10 9 m/sec and 

Mo = 4r x 10“ 7 henry /m, e 0 = Mo 1<: o 2 = (36jr) _l x 10 _n sec 2 /(henry.m) (5) 


t Other names are often used, see Table 1. 
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The condition /x « hq holds well for most practical purposes in such media as air and 
copper; in fact /x ai > = 1. 0000004 /x 0 and n copper = .99999/ip. 

The electrical field strength E is further related to the current density j by Ohm’s law: 

1 = »E (6) 

where g is the conductivity of the material. Again for an non-isotropic material g is 
generally a tensor which may also contain real and imaginary components; in which case 
the above relation becomes the generalized Ohm’s law. For good conductors g >> c; for 
bad conductors g « e. In free space, g = 0. 

8.8 Maxwell Equations in Terms of E and B 

To pass to the four-potential considered in Section 4 it is convenient to express Maxwell’s 
equations in terms of the electrical field strength E and the magnetic flux B. In fact this 
is the pair most frequently used in electromagnetic work that involve arbitrary media. On 
eliminating D and H through the constitutive equations (3) we obtain 


B + V x E = 0, 

VxB— fie E = /ij, 

II 

> 

V • B = 0. 


(7) 


The second equation assumes that e is independent of time; otherwise cE = c dE/dt should 
be replaced by d(eE)/dt. In charge-free vacuum the equations reduce to 


1 • 

B + V x E = 0, VxB = 0, 

c o 

V • E = 0, V • B = 0. 


( 8 ) 


8.4 The Electromagnetic Potentials 


The electric scalar potential $ and the magnetic vector potential A are introduced by the 
definitions 

(9) 


E = — — A, B = V x A. 


This definition satisfies the two homogeneous Maxwell equations in (7). The definition of 
A leaves its divergence V • A arbitrary. We shall use the Lorentz gauge [13] 


V • A + /xc$ = 0. 


( 10 ) 


With this choice the two non- homogeneous Maxwell equations in terms of $ and A separate 
into the wave equations 


V 2 $ — /xc$ = — p/e, V 2 A — fieA = — /xj. 


( 11 ) 
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4. THE ELECTROMAGNETIC FOUR-POTENTIAL 


Maxwell’s equations can be presented in a compact manner* in the four-dimensional space- 
time defined by the coordinates 


x\ = x, xa = y, 13 = 2 , x 4 = tcf 


( 12 ) 


where Xi, x 3 , x 3 are spatial Cartesian coordinates, i 7 = —1 is the imaginary unit, and 
e = 1 fyfixi is the speed of EM waves in the medium under consideration. In the sequel 
Roman subscripts will consistently go from 1 to 4 and the summation convention over 
repeated indices will be used unless otherwise stated. 

4-1 The Field Strength Tensor 

The unification can be expressed most conveniently in terms of the field-strength tensor 
F, which is a four-dimensional antisymmetric tensor constructed from the components of 
E and B as follows: 


0 

Fl2 

F 13 

Fu) 


° 

cB 3 

—cB 2 

~iE 1 > 

—Fia 

0 

Fas 

F 24 

d = 0 

— cB 3 

0 

cB 1 

-iEa 

—F 13 

—Fa 3 

0 

F 34 

cB a 

—cB\ 

0 

—iE z J 

-*14 

—Fax 

-^34 

0 J 



iEa 

iE z 

0 J 


(13) 


Here 0 is an adjustment factor to be determined later. Similarly, introduce the four-current 
vector J as 


J = 




w'l ' 


' Mi ] 

cfxja 

> = 0c < 

Ms \ 

Cfljz 

Ms I 

'Pi* > 


v iy/JTJep J 


(14) 


Then, for arbitrary 0, the non-homogeneous Maxwell equations, namely V x B — peE = //j 
and V • E = p/t, may be presented in the compact “continuity” formf 


dFa 

dx k 


(15) 


The other two Maxwell equations, V • B = 0 and V xE + B = 0 , can be presented as 


dF\k dFmi + dFkm 

dx m dxk dxi 


( 16 ) 


where the index triplet (i,j,k) takes on the values (1,2,3), (4,2,3), (4,3,1) and (4,1,2). 


* A form compatible with special relativity, 
t The covariant form of these two equations. 


4.2 The Four-Potential 

The EM “four-potential” 0 is a four-vector whose components are constructed with the 
electric and magnetic potential components of A and 


'0 r 

[ cA x 'j 

02 

def 1 cA>2 1 

03 

] CA3 j 

k 04 > 

l *• J 


4> = P< 


It may then be verified that F can be expressed as the four-curl of <f > , that is 

_ d<t> k d<f>i 

* “a ry I 

OXi OX k 

or in more detail and using commas to abbreviate partial derivatives: 


° 

01,2 — 02,1 

01,3 — 03,1 
\ 01,4 — 04,1 


02,1 — 01,2 

0 

02.3 — 03,2 

02.4 — 04,2 


03.1 — 01,3 

03.2 — 02,3 
0 

03,4 — 04,3 


04.1 — 01,4 ^ 

04.2 — 02,4 

04.3 ~ 03,4 

0 J 


(17) 


(18) 


(19) 


4.S The Lagrangian 

With these definitions, the basic Lagrangian of electromagnetism can be stated as$ 


L = \F ik F ik - Jib = {p* - Jib 

= iP 2 (c 2 b 2 - e J ) - ^-(j',^1 + / 2 a 2 +J3A3 - *>$), 


in which 

£ 2 = B T B = £ 2 + £ 2 + £ 2 , E 2 = E r E = El + El + E\. (21) 

Comparing the first term with the magnetic and electric energy densities [2,19,20] 

u m = ^B r H = -B 2 , Ue = iD T E=i £ E 2 , (22) 

Zfi 

we must have /? 2 c 2 = /3 2 /(pie) = l//x, from which 

= yfi. (23) 


f Lanczos [l 2] presents this Lagrangian for free space, but the expression (24) for an arbitrary 
material was found in none of the textbooks on electromagnetism listed in the References. 
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Consequently the required Lagrangian is 


L — — .B 2 “ — (ji^i + 32 A 2 + jzAs — pQ). 


(24) 


The associated variational form is 


R = 



LdVdt 


(25) 


where V is the integration volume considered in the analysis. In theory V extends over 
the whole space, but in the numerical simulation the integration is truncated at a distant 
boundary or special devices are used to treat the decay behavior at infinity. 


4-4 The Four- Field Equations 


On setting the variation of the functional (24) to zero we recover the field equations (15- 
16). Taking the divergence of both sides of (15) and observing that F is an antisymmetric 
tensor so that its divergence vanishes we get 


OJi 

— =cfi(V-j + p) =0, (26) 

The vanishing term in parenthesis is the equation of continuity, which expresses the law 
of conservation of charge. The Lorentz gauge condition (10) may be stated as V - <f> = 0. 
Finally, the potential wave equations (11) may be expressed in compact form as 


n<t>i = -j< (27) 

where □ denotes the “four-wave-operator” , also called the D’AIembertian: 

def a 2 a 2 a 2 a 2 a 2 

dikdxk dx\ + dx\ + dx\ c 2 dt 2 ' 

Hence each component of the four-potential </> satisfies an inhomogeneous wave equation. 
In free space, J,- = 0 and each component satisfies the homogeneous wave equation. 
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5. THE AXISYMMETRIC TEST EXAMPLE 


The simplest example for testing the finite element formulation based on the four-potential 
variational principle is provided by the axisymmetric magnetic field generated by a uniform, 
steady current flowing through a straight, infinitely long conducting wire of circular cross 
section. In the present Section we derive expressions for the magnetostatic fields outside 
and within the conductor. These analytical solutions will be later compared with the finite 
element numerical solutions. 

5.1 The Fret-Space Magnetic Field 

To take advantage of the axisymmetric geometry we choose a cylindrical coordinate system 
with the wire centerline as the longitudinal s-axis. The vector components in the cylindrical 
coordinate directions r, 9 and z are denoted by 

Aj, B i, E i in the r direction 
A 2 , i? 2 > E 2 in the 9 direction 
A 3 , B 3 , Ez in the z direction 

The electromagnetic fields will then vary in the radial direction (r) but not in the angular 
(0) and axial ( z ) directions. Similarly, the current density that flows in the wire has only 
one nonzero component acting in the positive or negative z direction; conventionally we 
select the positive direction. 

In Cartesian coordinates the radial component of the electrostatic potential in free space 
can be calculated from the expression (see, e.g., [2,10,18,19,20]) 

A - = A ^r«Lv\^ (29) 

where |r| is the distance between the elemental charge jz dV and the point in space at which 
we wish to find the field potential. The integral extends over the volume containing charges. 
This expression serves equally well in cylindrical coordinates. In fact, the transformation 
of z components will be one to one if the center of the systems coincide. 

As noted above the only non-vanishing component of the current vector is jz dS where dS 
is the elemental cross sectional area of the conductor and jz is the current density in the z 
direction. If dt represents the differential length of the wire, then f s jz dV = f s jz dS di = 
I dt = I dz and |r| = \/r 2 + z 2 . Substitution into Eq. (29) yields 



This integral diverges, but this difficulty can be overcome by taking the wire to have a 
finite length 2 L symmetric with respect to the field point, that is large with respect to its 
diameter. Integrating between — L and +L we get 
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Expanding this equation in powers of r/L and retaining only first-order terms gives 


A 3 




r + C. 


( 32 ) 


where C is an arbitrary constant. For subsequent developments it is convenient to select 
C = (noI/2ir) lnJ2r» where Rt is the “truncation radius” of the finite element mesh in the 
radial direction. Then 

/.. r\ / _ \ 

(33) 


*--(£) -(£) 


With this normalization A$ = 0 at r = Rt- Taking the curl of A gives the B field in 
cylindrical coordinates: 


f B1 l 

\B r ] 


B = V x A = { Bj > = < 

Be , 

► = < 

1*3 J 1 




1 dAa dA% 

~ ~fr 

dAi dAs 


> 

r o } 

► = < 

dA 3 \ 
I ( 

J 

l o J 


(34) 


1 f ” r^ 1 

It is seen that the only non-vanishing component of the magnetic flux density is 

dA 3 


B e = B 2 = h 0 H 2 = —gf = 


27rr " 


(35) 


This expression is called the law of Biot-Savart in the EM literature. 
5.2 Magnetic Field Within the Conductor 


Again restricting our consideration to the static case, we have from Maxwell’s equations 
in their integral flux form 


/ H ds= <[ n~ l B-da= f } - dS, (36) 

* C v C *5 

where C is a contour around the field point traversed counterclockwise with an oriented 
differential arclength da and dS is the oriented surface element inside the contour. The 
term for the electric field disappears in this analysis because E = 0. From before we know 
that the right hand side of Eq. (35) is equal to the normal component of the current that 
flows through the cross sectional area evaluated by the integral. In the free space case, this 
is the total current that flows through the conductor. But in the conductor the amount 
of current is a function of the distance r from the center. Again using I to represent the 
total current carried by the conductor, and R the radius of the conductor, and assuming 
an uniform current density j 3 = I/faR 2 ), the right hand side of (35) become_A 


Evaluating the left hand side of the integral and solving for B 2 gives: 


2irrn l B 2 = 


Bi = 


2jt R 2 ' 


Comparing with (34) we see that if fi = hq then B 2 is continuous at the wire surface r — R 
and has the value fiol l[2irR). But if p ^ no there is a jump [fi — no)I/[2irR) in B 2 . 

The magnetic potential A 3 within the conductor is easily computed by integrating —B 2 
with respect to r: 

ul r 2 

A ° = ~ ! kR> +c - < 39 > 

The value of C is determined by matching (33) at r = R, since the potential must be 
continuous. The result can be written 


[»" (* - 5*) ~ '“’k 


The preceding expressions (33)-(40) for A3 could also be derived in a somewhat more direct 
fashion by integrating the ordinary differential equation V 2 A 3 = r - 1 (< 9 (r< 9 A 3 /dr)dr) = 
fij3 to which the second of ( 11 ) reduces. 

6. FINITE ELEMENT DISCRETIZATION 

6.1 The Lagrangian in Cylindrical Coordinates 

To construct finite element approximations we need to express the Lagrangian (24) 


L =2^ B2 ~ i^ 2 -a r A-p*), (41) 

in terms of the potentials written in cylindrical coordinates. For B 2 we can use the 
expression of the curl (33) 


2 fldA 3 dA 2 \\ 

B IT ) + 


(dAi dAsV (1 

V dz dr ) \r 


d(rA 2 ) 1 dA\ 
dr r d6 


For E 2 we need the cylindrical-coordinate gradient formulas 

/ d$ , } \ 

(E 1 ) E r ) W + A ' 

E= £, = £, =- \1$* + aA, 

( Ez ) E ,J M 

1 "57 + As ' 


*•-*'*-(**2 


dA 2 y (d$ dAs 

+ ~dT j + [~d7 + ~dT 


(44) 



In the axisymmetric case, Ai — Az — 0 ; furthermore A x = A 3 is only a function of the 
radial distance from the wire. Therefore dAz/dO = dAzfdz = 0 . From symmetry consid- 
erations we also know that the electric field cannot vary in the 6 and z directions, which 
gives d$/dz = d$/d$ = 0 . Finally, the only nonvanishing current density component is 
js. Consequently the Lagrangian (41) simplifies to 



( 




- [jsAs - p«) . 


(45) 


6.2 Constructing EM Finite Elements 

To deal with this particular axisymmetric problem a two-node “line” finite element ex- 
tending in the radial r direction is sufficient. -In the following we deal with an individual 
element identified by superscript e. The two element end nodes are denoted by * and j. 
The electric potential # and the magnetic potential Az = A z are interpolated over each 
element as 

= A| = N^A5, (46) 

Here row vectors and contain the finite element shape functions for $ e and A$, 
respectively, which are only functions of the radial coordinate r: 


Ni = <**,(>•) JVi,.(r)>, N' a = ( N’ A ,(r) N' Aj (r) ) , (47) 

and column vectors and A 3 contain the nodal values of $ and Az, respectively, which 
are only functions of time t: 



(48) 


Substitution of these finite element assumptions into the Lagrangian (45) and then into 
Eq. (25) yields the variational integral as sum of element contributions R = where 


- -pN*$ e ) dV'dt. 


(49) 


where V e denotes the volume of the element. Taking the variation with respect to the 
element node values gives 


"* = /J' L < <a » t (Ir)" ^ a 3 + • a 3 - 73 (ni 


)' 


+ 


CL f- e (^r) r ^ + p (Ni 


(50) 


dV e dt. 
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On applying fixed-end initial conditions at t = to and t = t\ and the lemma of the calculus 
of variations, we proceed to equate each of the expressions in brackets to zero. From the 
first bracket we obtain for each element the following second-order dynamic equations for 
the magnetic potential at the nodes, which are purposedly written in a notation resembling 
the mass-stiffness-force equations of mechanics: 

M^Aj + K^AJ = r A , (51) 

where 

= = (52) 

tX = f h (N),) T rfK*. (53) 

Jv 

From the second bracket we obtain for the electric potential a simpler relation which does 
not involve time derivatives, ».e, is static in nature: 

K%9 e = fl, (54) 

where 



Assembling these equations in the usual way we obtain the semidiscrete master finite 
element equations: 

M a A 3 + K A A3 = f A , ^ 

= f«. j 


6.S The Static Case 

In time-independent problems, the term A3 disappears from (56) and the master finite 
element equations of electromagnetostatics become 


K a A 3 = f A > = f*. (57) 

If the current density and charge distributions are known a priori then these two equations 
may be solved separately. If only the charge distribution p is known then the second 
equation should be solved first to obtain the electric field E as gradient of the computed 
electric potential then the current density j can be obtained from Ohm’s law (6) and 
used to computed the force vector f A of the first equation, which is then solved for the 
magnetic potential. Conversely, if only the current density distribution is known a prior: 
the preceding steps are reversed. 

For the present test problem the current distribution is assumed to be known, and we shall 
be content with solving the first equation for the magnetic flux. 


14 



6.4 An Alternative Semidiscretization 

If upon setting the brackets of the variation (50) to zero we multiply them through by p 
and 1/e, respectively, the expressions for the mass, stiffness and force matrices become 


M - = L ? K - - L f; = L “^ TdV ’’ 



(58) 

The matrices M and K above are quite similar to the capacitance and reactance matrices, 
respectively, obtained in the potential analysis of acoustic fluids [7,8]. Another attractive 
feature of (58) is that = K« if the shape functions of both potentials coalesce, as is 
natural to assume. These advantages are, however, more than counterbalanced by the fact 
that “jump forces” contributions to fx and f« arise on material interfaces where p and 
e change abruptly, and the proper handling of such forces substantially complicates the 
programming logic. Note that this issue does not arise in the treatment of homogeneous 
acoustic fluids. 


6.5 Applying Boundary Conditions 

The finite element mesh is necessarily terminated at a finite size, which for the test prob- 
lem is defined as the truncation radius Rt alluded to in Section 5.1. In static calculations 
the material outside the FE mesh may be viewed as having zero permeability p, or, equiv- 
alently, infinite stiffness or zero potential. It follows that the potential value at the node 
located on the truncation radius may be prescribed to be zero. This is the only essential 
boundary condition necessary for this particular problem. 


7. NUMERICAL VALIDATION 
7.1 Finite Element Model 

The test problem consists of a wire conductor of radius R transporting a unit current 
density. For this problem the finite element mesh is completely defined if we specify the 
radial node coordinates r* = r* and r? = r* +1 for each element e. If the mesh contains N ec 
elements inside the conductor, those elements are numbered e = 1,2, ... N ec and nodes 
n = 1,2, ... N ee + 1 starting from the conductor center outwards. The first node (n = 1) 
is at the conductor center r = 0 and node n = N ee + 1 is placed at the conductor boundary 
r = R. The mesh is then continued with N e f elements into free space, with a double 
node at the counductor boundary. The last node is placed at r = Rt at which point the 
free space mesh is truncated; usually Rt = 4R to 512. Although the mesh appears to be 
one-dimensional, a typical element actually forms a “tube” of longitudinal axis z, internal 
radius r‘ and external radius r®, extending a unit distance along z. 
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Figure 1. Magnetic potential Aa vs. distance from center r, /x wir « = 10.0: finite 
element values (triangles) and analytical values (squares). 



Figure 2. Magnetic potential A3 vs. distance from center r, = 1.0: finite 

element values (triangles) and analytical values (squares). 
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0.060 1.000 2.000 T 3.000 4.000 5.000 

Figure S. Magnetic flux density Bi vs. distance from center r, p.,>, = 10.0: 
finite element values (triangles) and analytical values (squares). Values 
shown on the interface r = 1 with dark symbols have been extrapolated 
from element center values to display the jump more accurately; this 
extrapolation scheme has not been used elsewhere. 



Figure 4. Magnetic flux density Bi vs. distance from center r, = 1.0: finite 

element values (triangles) and analytical values (squares). 

17 





For the present study the magnetic potential was linearly interpolated in r, using the linear 
shape functions * 

Ni = <*(l-{) i(l + {)>, (59) 

where £ is the dimensionless isoparametric coordinate that varies from —1 at node * to 
+1 at node j. This interpolation provides for C° continuity of the potential inside the 
conductor and in free space. 

For the calculation of the element stiffnesses and force vectors, it was assumed that the 
permeability /x and the current density jz were uniform over the element. Then analytical 
integration over the element geometry gives 



where r m = ^ (r* + ry) is the mean radius and l = ry — r\ the radial length. For the test 
problem, /x is constant inside the conductor whereas outside it /x = /xo was assumed to 
be unity. The longitudinal current density is jz = //(ttR 2 ) inside the conductor whereas 
outside it jz vanishes. 


The master stiffness matrix and force vector were assembled following standard finite el- 
ement techniques. The only essential boundary condition was the setting of the nodal 
potential on the truncation boundary to zero, as explained in Section 6.5. The modified 
master equations were processed by a conventional symmetric skyline solver, which pro- 
vided the value of the magnetic potential at the mesh nodes. The magnetic flux density 
i? 2 ) which is constant over each element, was recovered in element by element fashion 
through the simple finite difference scheme 



This value is assigned to the center of element e. 


7.2 Numerical Results 


(61) 


The numerical results shown in Figures 1 through 6 pertain to a unit radius conductor 
( R = 1), with the external (free space) mesh truncated at Rt = 5. The element radial 
lengths r* — r\ were kept constant and equal to 0.25, which corresponds to 4 internal and 
16 external elements. 

The computed values of the potential Az are compared with the analytical solution given 
by Eqs. (33) and (40). As can be seen the agreement is excellent. The comparison between 
computed and analytical values of the magnetic flux density Bz shows excellent agreement 
except for the last element near the wire center, at which point the difference scheme 
(61) loses accuracy. The permeability of free space is conventionally selected to be unity. 
Figures 1, 3, and 5 illustrate the case where the wire permeability n w i re is set to 10.0, 
whereas Figures 2, 4, and 6 are for the case in which fx w ire is 1.0, that is, same as in 
free space. (The value of the susceptibility e does not appear in these magnetostatic 
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Figure S. Restriction of Figure 3 to r > R = 1 , /*««>• = 10.0, showing free space 
magnetic flux density in more detail. 



e.see i.64e 2.48e r 3.320 4.160 s.000 

Figure 6. Restriction of Figure A to r > R = 1, /*„,>, = 1.0, showing free space 
magnetic flux density in more detail. 
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computations.) Figures 1 and 2 show computed and analytical magnetic potentials. The 
slope discontinuity at r = 1 in Figure 1 is a consequence of the change in permeability p 
from the wire material to free space. Figures 3 and 4 show the computed and analytical 
magnetic flux densities. As discussed in Section 5.2, the jump at r = 1 in Figure 3 is due 
to the change in permeability p from the material to free space. Figures 5 and 6 show the 
computed and analytical magnetic flux densities in free space with more detail. Note that 
Figures 5 and 6 for r > 1 are identical; this is the expected result because, as shown in 
Section 5.1, the free-space magnetic flux field depends only upon the current enclosed by 
a surface integral around the wire and not on the details of the interior field distribution. 

In summary, the finite element model performed very accurately in the test problem and 
converged, as expected, to the analytical solution as the size of the elements decreased. 


8. CONCLUSIONS 

The results obtained in the one-dimensional steady-state case are encouraging, and appear 
to be extensible to two- and three-dimensional problems without major difficulties. The 
electric field remains effectively decoupled from the magnetic field except through Ohm’s 
law. Care must be taken, however, in modeling the forcing function terms so as to avoid 
the appearance of discontinuity-induced forces at physical interfaces. 

The next step in achieving the goal of a finite element model for a superconductor field is to 
study the time-dependent case, starting with harmonic currents and proceeding eventually 
to general transients. The code for this is currently written, but a suitable analytical 
solution for comparison with computed responses is still being developed. 

If encouraging results are obtained in the dynamic case, thermocoupling effects will be 
added to the code. References [3,17,22] discuss several different approaches applicable to 
various contexts ( e.g . eddy currents) and these will have to be investigated for suitability 
for capturing the couplings effects that are relevant to the superconducting problem. 

After modeling the coupling effects, the next step will be to model the superconducting 
fields. The feasibility of using the current model for superconductor applications is great, 
as the current density of a superconductor can be approximated by the standard current 
density multiplied by a constant squared. This constant is called the London penetration 
depth. Other analytical models that possess similar characteristics have been developed 
and are presented in Ref. [11]. 
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APPENDIX: COMPUTER PROGRAM 


This Appendix lists the computer program used to test the new electromagnetic elements on the 
axisymmetric test example. Sections of the program that pertain to the in-core skyline solver 
SKYFAC/SKYSOL and the command language reader TinyCllp are not listed here. Their source 
code is presented in the following publications: 

Felippa, C. A., Solution of Equations with Skyline-Stored Symmetric Coefficient Matrix, 
Computers & Structures, 5, 1975, pp. 13-25 

Felippa, C. A., A Command Reader for Interactive Programming, Engineering Computations, 
3, No. 3, 1985, pp. 203-238 


C-DECK AAVIRE 
OBLOCK FORTRAN 

program VIRE 
C 

Integer MUMEL. MUMNP, MDOF 

parameter (MUMEL- 100, MUMNP-MUMEL+1) 

parameter (MDOF-MUMNP) 
integer numel, numnp, ndof 

C 

character CCLVAL 

character status *80 

integer nodelm(2,MUMEL) , betag(MUMNP) 

double precision kmuelm(MUMEL) . kepselm(MUMEL) 
integer dlp(0:MDOF) 

double precision a(MUMNP) . b(MUMEL) 
double precision r (MUMNP) . f (MUMNP), f be (MUMNP) 
double precision sm(MD0F*3) 

double precision aex (MUMNP) , bex (MUMEL ) . f ex (MUMNP) 
double precision vl(MUMNP), v2(MUMEL) 
double precision kmu, keps, vrad, trad, inten 
integer nelwir, nelext 

C 

1000 continue 
C 

call MATERIAL (kmu. keps) 

call PRINTMAT (kmu. keps) 

call DIMENSIONS (wrad, trad) 

call PRINTDIM (vrad, trad) 

call CURRENT (inten) 

call PRINTCUR (inten) 

C 

1500 continue 
C 

call SUBDIVIDE (nelwir, nelext, numel, numnp, ndof) 

call PRINTSUB (nelwir, nelext, numel. numnp, ndof) 

call GENELEMS 
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$ (nalvir, nelext, kau, keps, nodala, kauala, kapaalo) 

call PRIHTELM (nuaal, nodala, kaualm, kapsala) 

call GEN NODES (nalvir, nalazt , vrad, trad, r) 

call GENBCTAG (nuanp, bctag) 

call PRINTIOD (nuanp. r. bctag) 

call GENEXACT (nuaal, nuanp, r, vrad, trad, 

$ kau, intan, aaz, bax) 

call ASSDBTF 

$ (nuaal, nodala, kauala, kapsala, 

$ nuanp, r. ndof , bctag, sa, dip, status) 

if (status .na. * *) go to 4000 

call SXYMUL (sa, ndof, dip, aax, fax, 0. vi, v2) 
call ASSEMFOR 

$ (nuaal, nodala, kauala, kapsala, intan, 

$ nuanp, r, vrad, trad, ndof, bctag, f, fbc, status) 

if (status .no. ' *) go to 4000 

call CLREAD (' Go ahaad and solva (j/n)? ', * *) 
if (CCLVAL(l) .na. *Y*) go to 4000 
call SKYCOV (sa, ndof, dip, fbc, a, status) 
if (status .na. * *) than 

call ERROR (*SEYC0V\ status) 
and if 

call PRINTSQL (nuanp, r, bctag, f, a, aax, fax) 
call MAGFIELD (nuaal, nodala, r, a, b) 

call PRINTMAG (nuaal, nodala, r, b, bax) 

C 

4000 continua 
C 

call CLREAD . (* Nav FE subdivision (y/n)? \ * *) 

if (CCLVAL(l) .aq. *Y*) go to 1500 
call CLREAD (‘ Nav problaa data (y/n)? ’, * *) 

if (CCLVAL(l) .aq. *Y*) go to 1000 
stop 
and 

OEND FORTRAN 
C=DECK ASSEMFOR 

C=*PURP0SE Assanbla forca vactor 
C-BLOCX FORTRAN 

subroutina ASSEMFOR 

$ (nuaal, nodala, kauala, kapsala, intan, 

$ nuanp. r. vrad, trad, ndof, bctag. f, fbc, status) 

intagar nuaal, nodaln(2,*), nuanp 

intagar ndof, bctag(ndof) 

intagar aldof(2) 

doubla pracision r(*) , kmuelo(*), kapsala(*) 
doubla pracision intan. vrad, trad 
doubla pracision f(*), fbc(*) 
character* (*) status 
doubla pracision ra(2) , f a(2) , au 
integer i, j, n, na 
C 

status *> * * 

do 1500 j - l.ndof 
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f(j) - 0.0 
1500 continue 
C 

do 3000 ne ■ 1 .numel 
do 2200 i - 1.2 
a ■ nodelm(i.ae) 
ro(i) - r(n) 
eldof (i) - a 
2200 continue 

bu ■ kmuala (ne) 

call FORCE (aa, r«, iataa, wrad, fa, status) 
if (atatua .aa. ' *} than 

call ERROR (‘ASSEMFOR’, atatua) 
and if 
C 

do 2500 i - 1.2 
j - eldof (i) 
f(j) - f(j) + fe(i) 

2500 continue 

C 

3000 continue 

do 4000 j - l.ndof 
fbc(j) - f(j) 

if (bctag(j) .na. 0) fbc(j) * 0.0 
4000 continue 
return 
and 

C=END FORTRAN 
C=DECX ASSEMSTF 

C=PURPOSE Aasambla master stiffness matrix 
C=BLOCX FORTRAN 

subroutine ASSS4STF 
$ (numel, nodelm, kmuelm, heps elm, 

$ numnp, r, ndof, betag, am, dip, status) 
character* (*) status 

integer numel, nodelm(2 , numel) , numnp 

integer ndof, bctag(ndof ) , dlp(0:ndof) 

double precision kmuelm (numel) , kaps elm (numel) 
double precision r(ndof), sm(*) 
double precision re(2), ame(2,2) 
integer eldof (2) 

integer i. j, k. ii, jj. n. ne 

C 

status * * * 

C 

call FORMDLP (numel, nodelm, ndof, betag, dip) 
do 2500 i - l.abs (dip (ndof)) 
sm(i) ■ 0.0 

2600 continue 
C 

do 4000 ne - 1, numel 
do 2200 i = 1 , 2 
n = nodelm(i.ne) 
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re(i) - r(a) 

•ldof(i) » a 
2200 continue 

' C 

cell STIFF (ne, re, kmuelm(ne), erne, status) 

If (status .m. * *) than 

cell ERROR ('ASSEMSTF*. status) 
end if 
C 

do 3600 i - 1.2 
ii - eldof(i) 
do 3600 J - i.2 
Jj - eldof(j) 
if (jj .1*. ii) then 

k - ebe(dlp(ii)) - ii ♦ jj 
sa(k) * sn(k) ♦ sne(i,j) 
end if 

3600 continue 

3600 continue 
C 

4000 continue 
C 

return 

end 

C-END FORTRAN 
C-DECK CURRENT 

C-PURPOSE Reed current intensity 
C-BLOCK FORTRAN 

subroutine CURRENT (inten) 
double precision DCLVAL 
double precision inten 

cell CLREAD (* Enter current intensity: ’, * ’) 

inten - DCLVAL(l) 

return 

end 

C-END FORTRAN 
C-DECK DIMENSIONS 

C-PURPOSE Reed problem dimensions (wire and truncation radius) 
C-BLOCK FORTRAN 

subroutine DIMENSIONS (wred, trad) 
double precision DCLVAL 
double precision wred, trad 

cell CLREAD (* Enter wire radius, trunc radius: *, * *) 

wred ■ DCLVAL (1) 

trad - DCLVAL (2) 

return 

end 

C-END FORTRAN 
C-DECK ERROR 

C-PURPOSE Fatal error termination subroutine 
C-BLOCK FORTRAN 

subroutine ERROR (name, message) 
character* (*) name, message 
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integer i , 1 
C 

1 “ len (mess age) 

do 1200 1 - len(message) ,1,-1 

if (bob sage (i :i) .ne. * ') go to 1300 
1 - i 

1200 continue 
1300 continue 

print *. * * 

print *, ’*** Fatal error condition detected ***' 

print *, messaged :1) 

print *, 'Error detected by name 

atop '*** Error stop ***' 

end 

C=END FORTRAN 
CODECS FORCE 

OPURPOSE Compute node forces for axisymm EM element due to j 
C=BL0CK FORTRAN 

subroutine FORCE 

$ (ne, re, inten, wrad, fe, status) 

integer ne 

double precision re(2), inten. vrad, fe(2) 
character* (*) status 
double precision ri, r j , rm, fn 
C 

status - ' ' 

ri ■ red) 
rj - re (2) 

if (rj .le. ri) then 

write (status , * (A, 15) ') 

$ ‘FORCE: Negative or zero length, element', ne 
return 
end if 

rm ■ 0.5*(ri+rj) 

if (rm . It . wrad) then 

fn ■ (inten/ (3 . 141B9*wrad**2) )* (rj-ri) 

fe(l) * fn* (ri+ri+r j )/6. 

fe(2) * fn*(ri+rj+rj)/6. 

else 

fe(l) - 0.0 

fe(2) « 0.0 

end if 
return 
end 

C=END FORTRAN 
C=DECK FORMDLP 

C=PURP0SE Form diagonal location pointer (DLP) array 
C=BL0CK FORTRAN 

subroutine FORMDLP 

$ (numel, nodelm, ndof, betag, dip) 

integer numel, nodelm (2, numel ) , ndof 

integer bctag(ndof ) , dlp(0:*) 

integer i. j, k, n. ne, eldof(2), nsky 
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c 

do 1200 i ■ O.ndof 
dlp(i) - 0 

1200 continue 

do 2000 ne ■ l.nunel 
do 1600 i - 1.2 
& ■ nodolad.no) 
eldof(i) ■ n 
1600 continue 

do 1800 i - 1.2 
k - eldaf (i) 

do 1800 J - 1.2 

if ( oldof (J) .lo. k) then 
dlp(k) - aax(dlp(k) ,k-eldof (j)+l) 
end if 

1800 continue 

2000 continue 

do 2200 i - l.ndof 

dlp(i) ■ dlp(i-l) ♦ dlp(i) 

2200 continue 

nsky « aba (dip (ndof ) ) 

C 

print *(/** No of equations: * *,110) * .ndof 

print * ( * * Average bandwidth : ' * ,F12 . 1 ) * . float (nsky) /ndof 

print *(** Entries to store skyline: ** ,110) * .nsky 

print *(*• **)* 

C 

do 3000 i * l.ndof 

if (bctag(i) .ne. 0) dlp(i) » -abs(dlp(i)) 

3000 continue 
return 
end 

OEND FORTRAN 
ODECK GENBCTAG 

C=PURP0SE Generate potential BC data by fixing externmost node 
OBLOCX FORTRAN 

subroutine GENBCTAG (nunnp, betag) 
integer nunnp . bctag(*) 

integer n 

do 2000 n * l.nxwnp 
bctag(n) ■ 0 
2000 continue 

bctag(nunnp) ■ 1 

return 

end 

C~END FORTRAN 
C=DECX GENELEMS 

C=PURP0SE Generate eleaent data 
C=«BL0CX FORTRAN 

subroutine GENELB4S 

$ (nelvir, nelext, kmu, keps, nodeln. knueln. kepseln) 

integer nelvir, nelext 

integer nx, n. ne, nodela(2,*) 
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double precision kmu, keps, kmuelm(*). kepselm(*) 

n ■ 0 

ne * 0 

do 2000 nx ■ 1 .nelvir 
n « n + 1 
ne ■ ne ♦ 1 
nodeln(l ,ne) ■ n 
nodeln(2.ne) ■ n+1 
kmuelm(ne) ■ kmu 
keps elm (ne) - keps 
2000 continue 

do 3000 nx ■ 1 .nelext 
n ■ n+1 

ne a ne + 1 
nodelm(l.ne) - n 
nodeln(2.ne) ■ n+1 

kmuela(ne) ■ 1.0 

kepseln(ne) - 1.0 

3000 continue 
return 
end 

OEND FORTRAN 
ODECX GENEXACT 

C=PURPQSE Generate exact magnetic potential/! i eld solutions 
C=BL0CK FORTRAN 

subroutine GENEXACT (numel, numnp. r, vrrad, trad, 

$ kmu, inten, aex, bex) 

integer numel . nunnp 

double precision r(*). wrad, trad, kmu. inten 
double precision aex(numnp) , bex(numel) 
integer n, ne 
double precision c, rm 
C 

c » -(inten/(2*3.14169))*log(vrad/trad) 
do 2000 n = 1 , n uan p 

if (r(n) .It. vrad) then 

aex(n) = (kmu*inten/(4*3.141696))*(l .-(r(n)/wrad)**2) + c 

else 

aex(n) = -(inten/(2*3. 14169))*log(r(n)/trad) 
end if 

2000 continue 

do 3000 ne = 1 .numel 

rm * 0.6*(r(ne)+r(ne+l)) 

if (rm .le. vrad) then 



bex(ne) = 

else 

(kmu* int en/ ( 2*3 . 1 4 1 69) ) * (rm/wrad* * 2 ) 


bex(ne) = 
end if 

(inten/(2*3 . 14160) )/rm 

3000 

continue 



return 

end 


C=END 

FORTRAN 


C=DECK GEN NODES 
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OPURPOSE Generate noda data 
OBLOCK FORTRAN 

aubroutina GENIODES (nalvir, nelext , vrad, trad, r) 

integer nalvir, nalaxt 

lntagar n, a • 

doubla pracialoa vrad, trad, r(*) 

C 

n ■ 0 

do 2000 na ■ 1 , nalvir 
n ■ n ♦ 1 

r(n) « (ne-l)*vrad/nelvir 
2000 continue 

r(n+l) * vrad 
do 3000 na ■ 1 .nalaxt 
n ■ n ♦ 1 

r (n) ■ vrad ♦ (na-l)*(trad-vrad)/nalaxt 
3000 contlnua 

r(n+l) ■ trad 

raturn 

and 

C-END FORTRAN 
ODECX MAGFIELD 

OPURPOSE Computed nagnatie 11 aid (B) at alaaant cantar 

OBLOCX FORTRAN 

aubroutina MAGFIELD (nunal, nodelm, r. a, b) 
lntagar nunal, nodaln(2,nunal) 

doubla praclaion r(*). a(*) . b(nunal) 
lntagar na, nl, nj 

C 

do 2000 na ■ 1, nunal 
nl * nodaln(l.na) 

nj “ nodaln(2,na) 

b(ne) - - (a(n j ) -a(ni ) ) / (r(nj ) -r (ni) ) 

2000 contlnua 
raturn 
and 

C=£ND FORTRAN 
C=DECK MATERIAL 

OPURPOSE Raad material propartlaa 
OBLOCX FORTRAN 

aubroutina MATERIAL (kmu, kepa) 
doubla praciaion knu, kapa 

doubla praciaion DCLVAL 

call CLREAD (* Enter knu, keps lor wire: *, * *) 

knu - DCLVAL ( 1 ) 

kepa - DCLVAL (2) 

raturn 

end 

C=END FORTRAN 
C=OECK PRINTCUR 

C=PURP0SE Print current intenaity 
OBLOCX FORTRAN 

aubroutina PRINTCUR (inten) 
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double precision int«n 

print *(** Current intensity: ** ,F10. 3) '.inten 

return 

end 

OEND FORTRAN 
ODECK PRINTDIM 

OPURPOSE Print problea dimensions (vire and truncation radiua) 

OBLOCK FORTRAN 

subroutine PRINTDIM (vrad, trad) 

double precision vrad, trad 

print *(•* Vire radius: ** .F10. 3) • .vrad 

print *(“ Truncation radius : * * .F10.3) * .trad 

return 

end 

OEND FORTRAN 
ODECK PRINTELM 
OPURPOSE Print element data 
OBLOCK FORTRAN 

subroutine PRINTQJ4 (numel, nodelm, kmuelm, leaps elm) 
integer i. n, numel, nodelm(2,*) 
double precision kmuelm (*) , kepsalm(*) 

print *. * * 

print *. ‘Element Data* 

print *, * * 

print *, 

$ * Elam I J kmu keps* 

do 2000 n ■ 1, numel 

print ‘(3I6.2F0.3) * . n. (nodelm(i ,n) ,i«l .2) ,kmuelm(n) .kepselm(n) 
2000 continue 
return 
and 

OEND FORTRAN 
ODECK PRINTMAG 

OPURPOSE Print computed and exact magnetic field (B) 

OBLOCK FORTRAN 

subroutine PRINTMAG (numel, nodelm, r, b, bex) 
integer numel, nodelm (2, numel) 

double precision r(*) ,b(numel) .bex(numel) 
integer ne . ni. nj 

C 

print *, * • 

print *. ‘Magnetic Field* 

print *, * * 

print * , 

$ * El am r-centar Comp-B2 Exact -B2* 

do 2000 ne = 1 .numel 
ni = nodelm (l,ne) 
nj ■ nodelm(2,ne) 

print ‘ (15 , F10. 3, 2F1 1.4)‘,ne,0.5* (r(ni) -r(nj) ) . 

$ b(ne) .bex(na) 

2000 continue 
return 
end 
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C-END FORTRAN 
C-DECK PRINTMAT 

C-PURPOSE Print material properties uasd in problem 
C-BLOCK FORTRAN 

subroutine PRINTMAT (kmu, kept) 
double precision kmu, keps 

print *('* Rel. permeability of wire (vacuum-1) : • * .F10.3) * ,kmu 
print *(** Rel. permittivity of wire (vacuum*!) :** ,F1 0.3) *, keps 
return 
end 

C-END FORTRAN 
C-DECK PRINTNOD 
C-PURPOSE Print element data 
C-BLOCK FORTRAN 

subroutine PRINTNOD (numnp , r, bctag) 
integer n, mump, bctag(*) 
double precision r(*) 

print *, ’ * 

print *, * Node Data* 

print *, * 

print * . * Node r-coord bctag* 

do 2000 n - 1, numnp 

print *(I6,F10.3,I6) * , n,r(n) .bctag (n) 

2000 continue 
return 
end 

C-END FORTRAN 
C-DECK PRINTSOL 

C-PURPOSE Print computed and exact solution 
C-BLOCK FORTRAN 

subroutine PRINTSOL (numnp, r, bctag, f. a, aex, fex) 
integer n, numnp, bctag (numnp) 
double precision r (numnp) , f (numnp) 

double precision a(numnp) . aex(numnp) , f ax(numnp) 

print * , * * 

print *, ‘Computed Solution’ 

print *, * * 

print *, 

$ * Node r bctag Comp-for*, 

$ * Comp- A3 Exact -A3 Exact -for* 

do 2000 n - 1, numnp 

print *(I6,F10.3,I6,4F11 .4) ’ ,n. r(n) .bctag (n) . 

$ f (n) ,a(n) ,aex(n) .fex(n) 

2000 continue 
return 
end 

C-END FORTRAN 
C-DECK PRINT SUB 

C— PURPOSE Print subdivision data 
C-BLOCK FORTRAN 

subroutine PRINTSUB (nelwir, nelext, numel, numnp, ndof) 

integer nelwir, nelext, numel, numnp, ndof 

print *(** Subdivisions in wire :**,I6)*, nelwir 
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o o o 


Subdivision* in free apaca 
Number of alanant* 

Number of noda point* 
Numbar of dof* 


'MS)’, nalaxt 
' ’,10)*. numal 
**,I6)*, nuanp 
*\I6)\ ndof 


print •(** 
print ■(" 
print * ( " 
print •(** 
raturn 
and 

OERD FORTRAN 
ODECK SKYCOV 

OPURPOSE Cover routina for atiffna** solvar 
C- AUTHOR C. A. Falippa. March 1072 
OVEtSION November 1082 (Fortran 77) 

OTHISVERSXON Condanaad on November 86 for ME503 HV 
C*EQUIPMENT Machina indapandant 
OKEYVORDS solve skyline stiffnasa aquation 
C=BL0CK ABSTRACT 
C 

SKYCOV ia a cover routina that solves tha master 
atiffnaaa aquations 

K u ■ f 

SKYCOV calls SKYFAC to factor tha skyline- stored 
master stiffness matrix K. If tha factorization is 
successful SKYCOV than calls SKYSOL to solva for u. 


C 
C 
C 
C 
C 
C 
C 

C=END ABSTRACT 
C»BL0CK USAGE 
C 

Tha calling sequence is 


CALL SKYCOV (S. N, DLP, F. U, STATUS) 
Input arguments : 


N 

P 

DLP 


C 
C 

c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 

C=END USAGE 
C=BL0CK FORTRAN 
subroutine 

$ 


Skyline stored stiffness matrix 
Overwritten by factorization. 
Number of aquations 
Nodal forca vector 
Skyline diagonal location pointer 


Output arguments : 


U 

STATUS 


Computed displacements if no error detected. 
Status character variable, 
blank no error detected 
nonblank explanatory error message 


SKYCOV 
(a, n, dip, 


f, u, status) 


ARGUMENTS 


integer 


n, dlp(0: *) 
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non non 


double precision e(*). f(*), u(*) 
charset er* (*) status 

TYPE t DIMENSION 

integer idetez, negeig, if ail, NMAX 

persaeter (NMAX*3000) 

double precision aux(NMAX) , detef , delta, DOTPRD 
external DOTPRD 

LOGIC 


status - ‘ * 

if (n .gt. NMAX) then 

write (status, * (A. Id) *) ‘No. of equations exceeds * .NMAX 
return 
end if 

call SKYFAC 

$ (s, 0, n. n, dip, aux, DOTPRD, .true., .false.. 

$ 0, 0, 0.0, datef, idetex, negeig, ifail) 

if (ifail .gt. O) than 

vrita (status, * (A, Id, A) ') 

$ ‘Factorization aborted at equation ‘.ifail, 

$ * (natrix appears singular) ’ 

return 
end if 
C 

call SKYSOL 

$ (s, n, dip, DOTPRD, 0, 1. f. u. 0, 0, aux, delta) 

C 

return 

end 

OEND FORTRAN 
C-DECK STIFF 

O>PURP0SE Construct stiffness matrix of axisymmatrie EM element 
OBLOCX FORTRAN 

subroutine STIFF (ne, re, mu, s, status) 
integer ne 

double precision re(2), mu, s(2,2) 
character* (*) status 
double precision ri, r j , rl. rm 
status - ‘ ‘ 

ri ■ re(l) 

rj « re(2) 

rl » rj - ri 

if (rl .la. 0.0) then 

write (status, * (A, 15) *) 

$ ‘STIFF: Negative or zero length, element *,na 

return 
end if 

rm * O.B*(ri-^rj) 
s(i.l)= rm/(rl*eu) 

■ (2.2)» s(l.l) 
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■( 1 . 2 )- -«( 1 . 1 ) 

■(2.1)- a(l .2) 

return 

end 

OEND FORTRAN 
C-DECX SUBDIVIDE 
C-PURPOSE Read subdivision data 
OBLOCK FORTRAN 

subroutine SUBDIVIDE (nelwir, nelext, nuael, nuamp, ndof) 
integer nelvir, nelert. nuael, nuamp, ndof 
Integer ICLVAL 

call CLREAD (' Subdivisions in wire: *, * *) 
nelwir « ICLVAL (1) 

call CLREAD ( a Subdivisions in free space: *, * *) 

nelext - ICLVAL (1) 

nuaiel - nelwir + nelext 

numnp - nuael ♦ 1 

ndof - nuanp 

return 

end 

C-END FORTRAN 
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